//==============================================================================
// LSGAffineAlignment.cc
// 
// Linear Space Affine gap penalty Global Alignmnet. 
// 
//==============================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8
//

#include "SequenceAlignment.hh"

#define STARTX  0x00000001
#define ENDX    0x00000010
#define STARTY  0x00000100
#define ENDY    0x00001000

//------------------------------------------------------
// * GetAlignment
//------------------------------------------------------
// Recovers the best globally aligned subsequence of two 
// sequences.  
//
void 
LSGAffineAlignment::GetAlignment()
{
  // Allocate two columns for score table and two columns 
  // for intersect forwarding table ("c" table). 
  mScoreTable.Allocate(12,mRows-1); 
  Align(0,0,0,mRows-1,mCols-1);
}

//------------------------------------------------------
// * Align
//------------------------------------------------------
// Since the gap penalty is affine, there are special cases
// to handle when the top left corner or lower right corner
// are known to be the beginning of a run of gaps.  Thus
// scaseTL and scaseLR
//
//
void 
LSGAffineAlignment::Align(char scaseTL,char scaseLR, int x1,int y1,int x2,int y2)
{
  int w = x2-x1;
  int h = y2-y1;
  int u,v;
  char dir=scase;

  // Sanity check
  if ((w<0) || (h<0)){
	 cout<<"Invalid range.\n";
	 exit(1);
  }

  if ((w < 2) && (h < 2)){
	 HandleLimitingCase(x1,y1,x2,y2);
  }else{
	 u = ((x2-x1)/2)+x1;
	 v = AffineFindInterceptWithBestPath2(scaseTL,scaseLR,scaseUV,u,x1,y1,x2,y2);
	 Align(scaseTL,scaseUV,x1,y1,u,v);
	 Align(scaseUV,scaseLR,u,v,x2,y2);
  }
}

//------------------------------------------------------
// * InitScoreTableColumns
//------------------------------------------------------
//
void
LSGAffineAlignment::InitScoreTableColumns(int rows)
{
  mScoreTable[SPrevious][0] = 0;
  mScoreTable[IxPrevious][0] = 0;
  mScoreTable[IyPrevious][0] = 0;

  // Fill in initial column with gap scores. 
  int i;
  for(i = 1;i<rows;i++){
	 mScoreTable[SPrevious][i] = mGapScore+(i-1)*mGapExtensionScore;
  }

  // Fill in initial column with gap-extension scores. 
  for(i = 1;i<rows;i++){
	 mScoreTable[IxPrevious][i] = mGapScore+(i-1)*mGapExtensionScore;
	 mScoreTable[IyPrevious][i] = mGapScore+(i-1)*mGapExtensionScore;
  }  
}

//------------------------------------------------------
// * InitColumnTops
//------------------------------------------------------
//
void
LSGAffineAlignment::InitColumnTops(int xIdx)
{
  // Fill in the top of each column with gap scores.
  mScoreTable[SCurrent][0] = mGapScore+mGapExtensionScore*(xIdx-1);
  mScoreTable[IxCurrent][0] = mGapScore+mGapExtensionScore*(xIdx-1);
  mScoreTable[IyCurrent][0] = mGapScore+mGapExtensionScore*(xIdx-1);

  // Top row of "c" table can only come from top,left corner.
  mScoreTable[CsCurrent][0] = 0+STable;
  mScoreTable[CxCurrent][0] = 0+XTable;
  mScoreTable[CyCurrent][0] = 0+YTable;
}


//------------------------------------------------------
// * FillWithBestScore
//------------------------------------------------------
//
void
LSGAffineAlignment::FillWithBestScore(int xIdx,int yIdx,int &FromWhere)
{
  int diagScore,diagIx,diagIy,topIy,leftIx,leftScore,topScore;
  int diagMax,topMax,leftMax,Max;

  // -1 adjusts for offset between string and table coordinates
  int  matchScore = Score((char)mS[x1+xIdx-1],mT[y1+yIdx-1]);

  diagScore = mScoreTable[SPrevious][yIdx-1] + matchScore;
  diagIx = mScoreTable[IxPrevious][yIdx-1] + matchScore;
  diagIy = mScoreTable[IyPrevious][yIdx-1] + matchScore;

  topScore = mScoreTable[SCurrent][yIdx-1]  + mGapScore;
  topIy = mScoreTable[IyCurrent][yIdx-1]  + mGapExtensionScore;
  
  leftScore =  mScoreTable[SPrevious][yIdx]  + mGapScore;	
  leftIx =  mScoreTable[IxPrevious][yIdx]  + mGapExtensionScore;	
  
  diagMax=max(diagScore,max(diagIx,diagIy));
  
  topMax = max(topScore,topIy);

  leftMax = max(leftScore,leftIx);
  
  int AllMax = max(leftMax,max(diagMax,topMax));

  mScoreTable[SCurrent][yIdx]= diagMax;
  mScoreTable[IyCurrent][yIdx]=topMax;
  mScoreTable[IxCurrent][yIdx]=leftMax;

}

//------------------------------------------------------
// * InitCTables
//------------------------------------------------------
//
void
LSGAffineAlignment::InitCTables(int yIdx)
{
  mScoreTable[CsCurrent][yIdx] = yIdx+STable;
  mScoreTable[CxCurrent][yIdx] = yIdx+XTable;
  mScoreTable[CyCurrent][yIdx] = yIdx+YTable;
}

//------------------------------------------------------
// * FillCTables
//------------------------------------------------------
//
void
LSGAffineAlignment::FillCTables(int yIdx,int CsOrigin,CxOrigin,CyOrigin)
{
  // Fill in Cs table 
  if (CsOrigin == SCORE){
	 mScoreTable[CsCurrent][yIdx] = 
		mScoreTable[CsPrevious][yIdx-1];
  }else{
	 if (CsOrigin == IX){
		mScoreTable[CsCurrent][yIdx] = 
		  mScoreTable[CxPrevious][yIdx-1];
	 }else{
		mScoreTable[CsCurrent][yIdx] = 
		  mScoreTable[CyPrevious][yIdx-1];
	 }
  }
			  
  // Fill in Cx table 
  if (CxOrigin == SCORE){
	 mScoreTable[CxCurrent][yIdx] = 
		mScoreTable[CsPrevious][yIdx];
  }else{
	 mScoreTable[CxCurrent][yIdx] = 
		mScoreTable[CxPrevious][yIdx];
  }
  
  // Fill in Cy table 
  if (CyOrigin == SCORE){
	 mScoreTable[CyCurrent][yIdx] = 
		mScoreTable[CsCurrent][yIdx-1];
  }else{
	 mScoreTable[CyCurrent][yIdx] = 
		mScoreTable[CyCurrent][yIdx-1];
  }
}			

//------------------------------------------------------
// * SwapColumns
//------------------------------------------------------
//
void
SequenceAlignemnt::SwapColumns()
{
  Swap(SPrevious,SCurrent);
  Swap(IxPrevious,IxCurrent);
  Swap(IyPrevious,IyCurrent);
  Swap(CsPrevious,CsCurrent);
  Swap(CxPrevious,CxCurrent);
  Swap(CyPrevious,CyCurrent);
}


//------------------------------------------------------
// * GetTraceBack
//------------------------------------------------------
//
int
LSGAffineAlignment::GetTraceBack(int rows,
										  int CsOrigin,
										  int CxOrigin,
										  int CyOrigin,
										  int &dir)
{

  int SMax = mScoreTable[SPrevious][rows-1];
  int IxMax = mScoreTable[IxPrevious][rows-1];
  int IyMax = mScoreTable[IyPrevious][rows-1];
  int ScoreMax = max(SMax,max(IxMax,IyMax));
  
  int TraceBack;

  if (ScoreMax == SMax){
	 // NOTE: I return Cprevious because the last swap made 
	 // Cprevious the final column. 
	 TraceBack = mScoreTable[CsPrevious][rows-1];
  }else{
	 if (ScoreMax == IxMax){
		TraceBack = mScoreTable[CxPrevious][rows-1];
	 }else{
		TraceBack = mScoreTable[CyPrevious][rows-1];
	 }
  }

  // Subtract out the table identifier and return the direction 
  // character. 
  if (TraceBack >= YTable){
	 TraceBack = TraceBack - YTable;
	 dir = 'Y';
  }else{
	 if (TraceBack >= XTable){
		TraceBack = TraceBack - XTable;
		dir = 'X';
	 }else{
		TraceBack = TraceBack - STable; // Typically zero. 
		dir = 'S';
	 }
  }
}

//------------------------------------------------------
// * FindInterceptWithBestPath
//------------------------------------------------------
//
int
LSGAffineAlignment::FindInterceptWithBestPath(char scase,
															char BestScoreOrigin,
															int u,int x1,
															int y1,int x2,
															int y2)
{
	int rows = y2-y1+1;
	int cols = x2-x1+1;
	int i,j;
	char ScoreOrigin;

	// Adjust u to local sub-table coordinates 
	u = u - x1;

	// Fill columns with appropriate initial values
	InitScoreTableColumns(rows);

	// Compute optimum alignment score one column at a time for 
	// the given range.  
	int x,y,yIdx,xIdx = 0;
	do{
	  xIdx++;
	  
	  // Put initial values in top of columns, handling any 
	  // special cases for the beginning of this alignment
	  InitColumnTops(xIdx,scase);

	  for(yIdx = 1;yIdx < rows;yIdx++){

		 // Fill cell xIdx,yIdx with best score computed from neighbors
		 FillWithBestScore(xIdx,yIdx,ScoreOrigin);

		 // Update the "c" tables which will map back the intercept.
		 // There are three separate "c" tables, one for each of Score,Ix,Iy.
		 if (xIdx == u){
			InitCTables(yIdx);
		 }else{
			if (xIdx > u){
			  FillCTables(yIdx,ScoreOrigin);
			}
		 }
	  }
	  
	  // Switch role of current/previous columns 
	  SwapColumns();

	}while((xIdx + x1) < x2);  	

	// Find the origin of the best score, handling 
	// any special cases for the end of this alignment
	FindBestScoreOrigin(scase,BestScoreOrigin);
	
	TraceBack = GetTraceBack(rows,BestScoreOrigin);

	// ScoreTable contains the traceback position for the -local- 
	// sub-table.  Add y1 to map it back to the full table. 
	return(y1+TraceBack);
}
